discover and transform plastic waste makers index data

purpose of notebook

(-) test and play with advanced numerical EDA methods

todos: (-) …

information

name: makeovermonday_2021w22 link: https://data.world/makeovermonday/2021w22 title: 2021/W22: The Plastic Waste Makers Index Data Source: Minderoo from 2019

insights

1. correlation - most of the columns are highly correlated, that was to be expected, since most variables are depend on each other, e.g., rigid + flexible = total c production, total -> -rank rigid overall has a less strong correlation with the other variables, which might hint to there being a a different sub-population based on rigid production, flexible has a stronger correlation with total as rigid and total, since flexible is a far bigger contribution to total flexible has a stronger correlation with the overall production, than rigid, this is interesting and also might hint to rigid producers have a different market strategy than flexible producers
2. a Gini coefficient of 0.56 is quite high, which means that only a top few producers are responsible for a large amount of the total SUP waste contribution, in other words 10% are responsible for 44% of the waste

load packages

overview

head(plastic)

summary(plastic)

      rank        polymer_producer    no_of_assets   production_of_in_scope_polymers flexible_format_contribution_to_sup_waste rigid_format_contribution_to_sup_waste
 Min.   :  1.00   Length:100         Min.   : 0.00   Min.   : 0.200                  Min.   :0.000                             Min.   :0.000                         
 1st Qu.: 25.75   Class :character   1st Qu.: 3.00   1st Qu.: 0.500                  1st Qu.:0.100                             1st Qu.:0.100                         
 Median : 50.50   Mode  :character   Median : 6.00   Median : 0.900                  Median :0.200                             Median :0.200                         
 Mean   : 50.50                      Mean   :11.56   Mean   : 1.805                  Mean   :0.538                             Mean   :0.416                         
 3rd Qu.: 75.25                      3rd Qu.:12.25   3rd Qu.: 1.700                  3rd Qu.:0.500                             3rd Qu.:0.500                         
 Max.   :100.00                      Max.   :82.00   Max.   :11.600                  Max.   :4.700                             Max.   :4.500                         
 total_contribution_to_sup_waste total_waste_div_production
 Min.   :0.200                   Min.   :0.3000            
 1st Qu.:0.300                   1st Qu.:0.4300            
 Median :0.450                   Median :0.5000            
 Mean   :0.950                   Mean   :0.5834            
 3rd Qu.:0.925                   3rd Qu.:0.6900            
 Max.   :5.900                   Max.   :1.0000            

correlation
most of the columns are highly correlated, that was to be expected, since most variables are depend on each other, e.g., rigid + flexible = total c production, total -> -rank rigid overall has a less strong correlation with the other variables, which might hint to there being a a different sub-population based on rigid production, flexible has a stronger correlation with total as rigid and total, since flexible is a far bigger contribution to total flexible has a stronger correlation with the overall production, than rigid, this is interesting and also might hint to rigid producers have a different market strategy than flexible producers
%>% select(-rank, -total, -assets) can show a more clear picture by removing dependent variables
```r name = c(’’) df <- plastic %>% select(-polymer_producer, - total_waste_div_production) %>% mutate(rank = -rank) %>% # change sign of rank to make it increase with the dependent variables rename( assets = no_of_assets, product = production_of_in_scope_polymers, flexible = flexible_format_contribution_to_sup_waste, rigid = rigid_format_contribution_to_sup_waste, total = total_contribution_to_sup_waste)
library(corrplot) # correlation plots # https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html
cor <- cor(df) cor_mtest <- cor.mtest(df, conf.level = 0.99) # combining correlogram with significance test corrplot(cor, method = “number”, order = ‘hclust’, addrect = 3, p.mat = cor_mtest$p, insig = “pch”) # order = AOE, FPC, hclust + addrect ```
```r
corrplot(cor, p.mat = cor_mtest\(p, low = cor_mtest\)lowCI, upp = cor_mtest$uppCI, order = ‘hclust’, sig.level = 0.01, tl.pos = ‘d’, addrect = 3, rect.col = ‘navy’, plotC = ‘rect’, cl.pos = ‘n’, insig = “p-value”) ```

Lorenz curve & Gini coefficient

a Gini coefficient of 0.56 is quite high, which means that only a top few producers are responsible for a large amount of the total SUP waste contribution, in other words 10% are responsible for 44% of the waste

name = c('total_contribution_to_sup_waste', 'Producers')
df <- plastic %>% rename(value = total_contribution_to_sup_waste) %>% select(value) 

library(gglorenz) #transformations for plotting Lorenz curve, https://github.com/jjchern/gglorenz

lorenzcurve <- df %>% 
  ggplot(aes(value)) +
    stat_lorenz(desc = FALSE) +
    coord_fixed() +
    geom_abline(linetype = 'dashed') +
    theme_minimal() +
    hrbrthemes::scale_x_percent() +
    hrbrthemes::scale_y_percent() +
    # hrbrthemes::theme_ipsum_rc() +
    annotate_ineq(df$value) +
    ggtitle(paste("Lorenz curve for", name[1], sep=" ")) 
lorenzcurve <- ggplotly(lorenzcurve) %>% layout(yaxis = list(title = paste("cumulative percentage of", name[1], sep=" ")), xaxis = list(title = paste("cumulative percentage of", name[2], sep=" "))) 

lorenzcurve

name = c('flexible_format_contribution_to_sup_waste', 'rigid_format_contribution_to_sup_waste', 'Producers', 'flexible', 'rigid')
df <- plastic %>% rename(flexible = flexible_format_contribution_to_sup_waste, rigid = rigid_format_contribution_to_sup_waste) %>% select(flexible, rigid) %>% pivot_longer(cols = c(flexible,rigid))

library(gglorenz) #transformations for plotting Lorenz curve, https://github.com/jjchern/gglorenz

gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}

lorenzcurve <- df %>% 
  ggplot(aes(x = value, color = name)) +
    stat_lorenz(desc = FALSE) +
    coord_fixed() +
    geom_abline(linetype = 'dashed') +
    theme_minimal() +
    hrbrthemes::scale_x_percent() +
    hrbrthemes::scale_y_percent() +
    # hrbrthemes::theme_ipsum_rc() +
    annotate_ineq(filter(df, name == name[4])$value, y = 0.95, colour = gg_color_hue(2)[1]) +
    annotate_ineq(filter(df, name == name[5])$value, y = 0.90, colour = gg_color_hue(2)[2]) +
    ggtitle(paste("compare Lorenz curve of", name[1], "and", name[2], sep=" ")) 
lorenzcurve <- ggplotly(lorenzcurve) %>% layout(yaxis = list(title = paste("cumulative percentage of<br>", name[1], "<br>", name[2], sep="")), xaxis = list(title = paste("cumulative percentage of", name[2], sep=" "))) 

lorenzcurve

principal component analysis colored by self organizing map cluster

I need more knowledge how to work with and interpret SOM and PCA, maybe also not enough observations in data set https://iamciera.github.io/SOMexample/html/SOM_RNAseq_tutorial_part2a_SOM.html

name = c('polymer_producer')
df <- plastic %>% select(- total_waste_div_production, -rank, -no_of_assets, -total_contribution_to_sup_waste) %>% # removed variables which are depended on each other
  rename( product = production_of_in_scope_polymers, 
          flexible = flexible_format_contribution_to_sup_waste, 
          rigid = rigid_format_contribution_to_sup_waste)


library(kohonen) # functions to train self-organising maps (SOMs)

# setup for pca
scale_data <- as.matrix(t(scale(t(df[, !names(df) %in% name])))) # We need to normalize the data based on scale function because the variables are different scales; Normalization means subtracting mean from each observation and dividing with standard deviation. Check all the variables mean values are zero now
head(scale_data)

       product   flexible      rigid
[1,] 1.0838622 -0.1970659 -0.8867964
[2,] 1.0301899 -0.0633963 -0.9667936
[3,] 1.1172018 -0.3058262 -0.8113756
[4,] 0.6859477 -1.1474034  0.4614557
[5,] 1.1208971 -0.3202563 -0.8006408
[6,] 1.0995250 -0.2443389 -0.8551861

# principle component analysis
pca <- prcomp(scale_data, scale=TRUE)
summary(pca) 

Importance of components:
                          PC1     PC2       PC3
Standard deviation     1.6591 0.49736 4.375e-16
Proportion of Variance 0.9175 0.08246 0.000e+00
Cumulative Proportion  0.9175 1.00000 1.000e+00

# add back to original so everything is together
pca_scores <- data.frame(pca$x)
data_val <- cbind(df, pca_scores)

pca_plot <- ggplot(data_val, aes(x = PC1, y = PC2)) +
    geom_rug(alpha = 0.5) + # two 1d marginal distributions, display individual cases so are best used with smaller datasets
    geom_density_2d(alpha = 0.2, bins = 4) +# 2D kernel density estimation using MASS::kde2d() and display the results with contours
    geom_point(alpha = 0.75) + # point geom is used to create scatterplots
    theme_minimal()
pca_plot <- ggplotly(pca_plot) %>% layout()

pca_plot

# clustering is performed using the som() function on the scaled gene expression values.
set.seed(3)

# define a grid for the SOM and train
grid_size <- ncol(scale_data)
som_grid <- somgrid(xdim = grid_size, ydim = grid_size, topo = 'hexagonal')
som_model <- som(scale_data, grid = som_grid)
summary(som_model)

SOM of size 3x3 with a hexagonal topology and a bubble neighbourhood function.
The number of data layers is 1.
Distance measure(s) used: sumofsquares.
Training data included: 100 objects.
Mean distance to the closest unit in the map: 0.007.

# generate som plots after training
plot(som_model, type = 'mapping')

plot(som_model, type = 'codes')

# plot(som_model, type = 'counts')
# plot(som_model, type = 'dist.neighbours')
# plot(som_model, type = 'quality')
# plot(som_model, type = 'changes')

# further split the clusters into a smaller set of clusters using hierarchical clustering.
som_cluster <- cutree(hclust(dist(som_model$codes[[1]])), 2) # use hierarchical clustering to cluster the codebook vectors

plot(som_model, type="mapping", bgcol = som_cluster, main = "Clusters")
add.cluster.boundaries(som_model, som_cluster)

# attach the hierchal cluster to the larger dataset data_val.
gridSquare <- grid_size * grid_size
som_clusterKey <- data.frame(som_cluster)
som_clusterKey$unit_classif <- c(1:gridSquare)
data_val <- cbind(data_val,som_model$unit.classif,som_model$distances)  %>% rename(unit_classif = 'som_model$unit.classif', distances = 'som_model$distances')
data_val <- merge(data_val, som_clusterKey, by.x = "unit_classif" )
head(data_val)

# plot pca with colored clusters
pcasom_plot <- ggplot(data_val, aes(x = PC1, y = PC2, color = factor(som_cluster), text = polymer_producer)) +
    geom_rug(alpha = 0.5) + # two 1d marginal distributions, display individual cases so are best used with smaller datasets
    geom_point(alpha = 0.75) + # point geom is used to create scatterplots
    theme_minimal()
pcasom_plot <- ggplotly(pcasom_plot) %>% layout()

pcasom_plot

# two variables, continuous x, continuous y, show trend and distribution
name = c('production_of_in_scope_polymers', 'total_contribution_to_sup_waste')
df <- merge(plastic, data_val, by.x = 'polymer_producer')
df <- df %>% rename(x = production_of_in_scope_polymers, y = total_contribution_to_sup_waste, cluster = som_cluster, text = polymer_producer) %>% select(x, y, cluster, text) 

# https://ggplot2.tidyverse.org/reference/geom_smooth.html
point_plot <- df %>%
  ggplot(aes(x = x, y = y, color = factor(cluster))) +
    # geom_jitter(alpha = 0.5, size = 1) +
    geom_rug(alpha = 0.5) + # two 1d marginal distributions, display individual cases so are best used with smaller datasets
    geom_density_2d(alpha = 0.2, bins = 4) +# 2D kernel density estimation using MASS::kde2d() and display the results with contours
    geom_smooth(fill = "grey90") + # aids the eye in seeing patterns in the presence of overplotting
    geom_point(aes(text = text), alpha = 0.75) + # point geom is used to create scatterplots
    theme_minimal() +
    ggtitle(paste("trend of", name[2], "over", name[1], sep=" ")) 

Ignoring unknown aesthetics: text

point_plot <- ggplotly(point_plot) %>% layout(xaxis = list(showticklabels = FALSE))

`geom_smooth()` using method = 'loess' and formula 'y ~ x'

x_density_plot <- df %>%
  ggplot(aes(x = x, color = factor(cluster))) +
    stat_density(geom="line") + # draws kernel density estimate, which is a smoothed version of the histogram
    # geom_histogram(binwidth = 1) +
    theme_minimal() 
x_density_plot <- ggplotly(x_density_plot) %>% layout(yaxis = list(showticklabels = FALSE, showgrid = FALSE), xaxis = list(showticklabels = FALSE, showgrid = FALSE))

y_density_plot <- df %>%
  ggplot(aes(x = y, color = factor(cluster))) +
    stat_density(geom="line") + # draws kernel density estimate, which is a smoothed version of the histogram
    # geom_histogram(binwidth = 1) +
    coord_flip() +
    theme_minimal() 
y_density_plot <- ggplotly(y_density_plot) %>% layout(yaxis = list(showticklabels = FALSE, showgrid = FALSE), xaxis = list(showticklabels = FALSE, showgrid = FALSE))

# https://ggplot2.tidyverse.org/reference/geom_quantile.html
qualtile_plot <- df %>%
  ggplot(aes(x = x, y = y, color = factor(cluster))) +
    geom_quantile(alpha = 0.8) + # fits a quantile regression to the data and draws the fitted quantiles with lines
    theme_minimal() 
qualtile_plot <- ggplotly(qualtile_plot) %>% layout(yaxis = list(showticklabels = FALSE, showgrid = FALSE))

Smoothing formula not specified. Using: y ~ x
Smoothing formula not specified. Using: y ~ x

# merge figures into one plot, via subplots, https://plotly-r.com/arranging-views.html
sub1 <- subplot(x_density_plot, plotly_empty(), point_plot, y_density_plot, nrows = 2, margin = 0, heights = c(0.15, 0.85), widths = c(0.9, 0.1), shareX = TRUE, shareY = TRUE, titleX = FALSE, titleY = FALSE) %>% layout()

No trace type specified and no positional attributes specifiedNo trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plotly.com/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Can only have one: config

sub2 <- subplot(qualtile_plot, plotly_empty(), margin = 0, widths = c(0.9, 0.10), titleX = FALSE, titleY = FALSE) %>% layout()

No trace type specified and no positional attributes specifiedNo trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plotly.com/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Can only have one: config

fig <- subplot(sub1, sub2, nrows = 2, margin = 0, heights = c(0.8, 0.2), shareX = TRUE) %>% layout(xaxis = list(title = name[1]), yaxis = list(title = name[2]))

Can only have one: config

  
fig